rm(list=ls(all=TRUE))

library(foreign) #import data
library(MASS) # some regression tools
library(spdep) #to look for spatial autocorrelation
library(maptools) #spatial tools
library(spatstat) #even more spatial tools
library(sandwich) #cluster; robust se
library(lmtest) # same as above
library(rms) # logit models

# get the data
setwd("~/Desktop")
data <- read.csv("RAP_data.csv")

# base models (ols)
h1 <- lm(severity ~ dir_int + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data)
summary(h1)
h2a <- lm(severity ~ purs_reb + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data)
summary(h2a)
h2b <- lm(severity ~ deter_ext + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data)
summary(h2b)
h2c <- lm(severity ~ int_coerce + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data)
summary(h2c)
h3 <- lm(severity ~ spillover + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data)
summary(h3)
h4a <- lm(severity ~ dir_int + territor + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data)
summary(h4a)
h4b <- lm(severity ~ int_coerce + territor + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data)
summary(h4b)

## build spatial weights (with some variations)
coordinates(data) <- ~ longitude + latitude
coord.matrix <- coordinates(data)
IDs <- row.names(as(data, "data.frame"))
nnb.four <- knn2nb(knearneigh(coord.matrix, k=4), row.names=IDs)
nnb.eight <- knn2nb(knearneigh(coord.matrix, k=8), row.names=IDs)
nnb.twelve <- knn2nb(knearneigh(coord.matrix, k=12), row.names=IDs)

wgt.four <- nb2listw(nnb.four, style="W")
wgt.eight <- nb2listw(nnb.eight, style="W")
wgt.twelve <- nb2listw(nnb.twelve, style="W")

## spatial tests for dependence
lm.morantest(h1, wgt.four)
lm.morantest(h2a, wgt.four)
lm.morantest(h2b, wgt.four)
lm.morantest(h2c, wgt.four)
lm.morantest(h3, wgt.four)
lm.morantest(h4a, wgt.four)
lm.morantest(h4b, wgt.four)

## Ok, so what type of model?
lm.LMtests(h1, wgt.four, test="all") #Lag
lm.LMtests(h2a, wgt.four, test="all") #Lag
lm.LMtests(h2b, wgt.four, test="all") #Lag
lm.LMtests(h2c, wgt.four, test="all") #Lag
lm.LMtests(h3, wgt.four, test="all") #Lag
lm.LMtests(h4a, wgt.four, test="all") #Lag
lm.LMtests(h4b, wgt.four, test="all") #Lag

# do the lag models
h1.lag <- lagsarlm(severity ~ dir_int + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data, wgt.four)
summary(h1.lag)
h2a.lag <- lagsarlm(severity ~ purs_reb + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data, wgt.four)
summary(h2a.lag)
h2b.lag <- lagsarlm(severity ~ deter_ext + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data, wgt.four)
summary(h2b.lag)
h2c.lag <- lagsarlm(severity ~ int_coerce + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data, wgt.four)
summary(h2c.lag)
h3.lag <- lagsarlm(severity ~ spillover + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data, wgt.four)
summary(h3.lag)
h4a.lag <- lagsarlm(severity ~ dir_int + territor + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data, wgt.four)
summary(h4a.lag)
h4b.lag <- lagsarlm(severity ~ int_coerce + territor + log(capratio) + demdem + cwpceyrs + factor(StYear), data=data, wgt.four)    
summary(h4b.lag)

# print the results
require(stargazer)
stargazer(h1.lag, h2a.lag, h2b.lag, h2c.lag, h3.lag, h4a.lag, h4b.lag)

# do logit models (onset of interstate war)
dir.int <- lrm(midwar ~ dir_int + ln_cap + demdem + cwpceyrs, x=T, y=T, data=data)
dir.int <- robcov(dir.int, cluster=data$cdyad)
dir.int
reb <- lrm(midwar ~ purs_reb + ln_cap + demdem + cwpceyrs, x=T, y=T, data=data)
reb <- robcov(reb, cluster=data$cdyad)
reb
deter <- lrm(midwar ~ deter_ext + ln_cap + demdem + cwpceyrs, x=T, y=T, data=data)
deter <- robcov(deter, cluster=data$cdyad)
deter
spill <- lrm(midwar ~ spillover + ln_cap + demdem + cwpceyrs, x=T, y=T, data=data)
spill <- robcov(spill, cluster=data$cdyad)
spill
dir.int.ter <- lrm(midwar ~ dir_int + territor + ln_cap + demdem + cwpceyrs, x=T, y=T, data=data, se.fit=T)
dir.int.ter <- robcov(dir.int.ter, cluster=data$cdyad)
dir.int.ter
ic <- lrm(midwar ~ int_coerce + ln_cap + demdem + cwpceyrs, x=T, y=T, data=data, se.fit=T)
ic <- robcov(ic, cluster=data$cdyad)
ic
ic.ter <- lrm(midwar ~ int_coerce + territor + ln_cap + demdem + cwpceyrs, x=T, y=T, data=data, se.fit=T)
ic.ter <- robcov(ic.ter, cluster=data$cdyad)
ic.ter

# print results 
stargazer(dir.int, reb, deter, ic, spill, dir.int.ter, ic.ter)

# do predicted probabilities
base <- Predict(dir.int.ter, dir_int=0, territor=0, ln_cap=mean(data$ln_cap), demdem=0, cwpceyrs=mean(data$cwpceyrs), fun=plogis)
p1 <- Predict(dir.int.ter, dir_int=1, territor=0, ln_cap=mean(data$ln_cap), demdem=0, cwpceyrs=mean(data$cwpceyrs), fun=plogis)
p2 <- Predict(reb, purs_reb=1, ln_cap=mean(data$ln_cap), demdem=0, cwpceyrs=mean(data$cwpceyrs), fun=plogis)
p3 <- Predict(deter, deter_ext=1, ln_cap=mean(data$ln_cap), demdem=0, cwpceyrs=mean(data$cwpceyrs), fun=plogis)
p4 <- Predict(spill, spillover=1, ln_cap=mean(data$ln_cap), demdem=0, cwpceyrs=mean(data$cwpceyrs), fun=plogis)
p5 <- Predict(dir.int.ter, dir_int=0, territor=1, ln_cap=mean(data$ln_cap), demdem=0, cwpceyrs=mean(data$cwpceyrs), fun=plogis)

# build dataframe for plotting
names <- c("Baseline", "Direct Intervention", "Pursue Rebels", "Deter Externalization", "Spillover Event", "Territorial Dispute")
prob <- c(base$yhat, p1$yhat, p2$yhat, p3$yhat, p4$yhat, p5$yhat)
lower <- c(base$lower, p1$lower, p2$lower, p3$lower, p4$lower, p5$lower)
upper <- c(base$upper, p1$upper, p2$upper, p3$upper, p4$upper, p5$upper)

plot.df <- data.frame(name=names, prob=prob, lower=lower, upper=upper)

# make the plot
png(file="~/Desktop/RAP Revision/warprobs.png", width=5000, height=3296, res=380)
plot(plot.df$prob, type="n", ylim=c(-0.01,0.15), ylab="Probability of War Onset", xlab="Type of Dispute", xaxt="n")
axis(1, at=1:6, labels=c("basline", "direct intervention", "pursue rebels", "deter externalization", "spillover", "territorial dispute"))
abline(h=-0.01, lty=2, col="gray")
abline(h=0, lty=2, col="gray")
abline(h=0.01, lty=2, col="gray")
abline(h=0.02, lty=2, col="gray")
abline(h=0.03, lty=2, col="gray")
abline(h=0.04, lty=2, col="gray")
abline(h=0.05, lty=2, col="gray")
abline(h=0.06, lty=2, col="gray")
abline(h=0.07, lty=2, col="gray")
abline(h=0.08, lty=2, col="gray")
abline(h=0.09, lty=2, col="gray")
abline(h=0.10, lty=2, col="gray")
abline(h=0.11, lty=2, col="gray")
abline(h=0.12, lty=2, col="gray")
abline(h=0.13, lty=2, col="gray")
abline(h=0.14, lty=2, col="gray")
abline(h=0.15, lty=2, col="gray")
lines(plot.df$prob, col="red", lty=1, lwd=2)
lines(plot.df$upper, col="gray23", lty=2, lwd=2)
lines(plot.df$lower, col="gray23", lty=2, lwd=2)
points(plot.df$prob, pch=21, bg="tomato2", col="black", cex=1.5)
points(plot.df$upper, pch=21, bg="black", col="black", cex=0.5)
points(plot.df$lower, pch=21, bg="black", col="black", cex=0.5)
dev.off()
